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ABSTRACT 


This report describes a multi phase self-adaptive predictor 
corrector type algorithm. This type of algorithm enables the solu- 
tion of highly nonlinear structural responses including kinematic, 
kinetic and material effects as well as potential pre/postbuckl i ng 
behavior. The hierarchy of the strategy is such that three main 
phases are involved. The first features the use of a warpable 
hyperel 1 i pti c constraint surface which serves to upperbound depen- 
dent iterate excursions during successive Incremental Newton 
Raphson type iterations. The second corrector phase uses an energy 
constraint to scale the generation of successive iterates so as to 
maintain the appropriate form of local convergence behavior. The 
third involves the use of quality of convergence checks which 
enable various self-adaptive modifications of the algorithmic 
structure when necessary. Such restructuring is achieved by tight- 
ening various conditioning parameters as well as switch to different 
algorithmic levels so as to improve the convergence process. 

Several numerical experiments illustrate the capabilities of the 
procedure to handle varying types of nonlinear structural behavior. 
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1 . INTRODUCTION 

One of the central features in the development of finite 
element computer programs for nonlinear analysis is the proper 
selection of solution algorithms. The nature of structural 
nonlinearities is generally quite diverse when both kinematic and 
material effects are included. Specifically, for static problems, 
such effects give rise to nonlinear algebraic equations which may 
possess path dependent multiple solutions. In this context, 
the quest for reliable as well as computationally efficient sol- 
utions to such problems is a very demanding task. 

Solution procedures for nonlinear problems have been discussed 

by a multitude of authors this direction, the mature 

works of Bergan et al . ^ Ri ks^^*^^and Cri sf i el d give a good 

overview of much of the progress made to date. As can be seen from 

rg-12 T 

these papers , unlike linear problems, it is extremely 

difficult to develop a single methodology of general validity 
which can be used to handle the diversity of potential structural 
problems. Since the formulation of the problem and hence the 
associated computer coding architecture is strongly dependent on 
the algorithmic approach taken, generally most general purpose 
(GP) nonlinear finite element (FE) codes have adopted one part- 
icular methodology through which the nonlinear problem is 
sol ved^^^“^*^^ . In this context, generally some variant of the 
Incremental Newton Raphson (INR) approach has been chosen^^^“^^ ^ . 
While the INR procedure is perhaps the most powerful of the 
iterative solution techniques, it is subject to several short- 
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comings. The more important of these can be categorized as follows 

1) Inefficiencies associated with update requirements; and 

2) Sensitivities/anomalous convergence characteristics in 
the neighborhood of turning points (zones of changing 
curvature definiteness), bifurcations, "shallow" curvature, 
snap through behavior, etc. 

r X B — 1 ft 1 

While the recently advocated pseudo update procedures 
provide a partial answer to the computational inefficiencies 
associated with updating, no real improvement is achieved in 
category 2) problems nor is it clear what happens in path dependent 
and/or multiple solution situations. 

To overcome the sensitivities associated with the use of the 
INR algorithm in the vicinity of turning points several approaches 
have been advocated, in particular: 

a) Use of deflection control^^^^; 

b) Rotation of solution space via introduction of auxiliary 
stiffness^^^^ ; 

c) Switch from step- i terati ve to pure Euler-Cauchy type 

f9 1 

incrementations initiated via curvature monitoring ; and 

d) Use of constraints to control successive dependent iterate 
excursi ons^^*^”^^^ . 

Since the main sensitivities/anomalous behavior of the INR type 
algorithms appears to be the generation of excessive iterate 
excursions in neighborhoods of turning point, shallow curvature 
etc., the constrained approach advocated by d) ^^0-12] gppg-^rs to 
be the best choice for use in general purpose (GP) codes. 
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The difficulty of the foregoing approaches lies in the fact 
that there is no automatic correction features associated with the 
algorithms wherein, as the solution proceeds, its quality^^^^ is 
monitored so as to enable the appropriate automatic corrective 
action to be taken. In this context, this report develops a self- 
adaptive type predictor-corrector algorithmic strategy. The hier- 
archy of the strategy is such that the predictor phase consists 
of a constrained type INR algorithm (CINR) which is employed 
to tunnel into the solution space. It features the use of a 
warpable hyperel 1 i pti c constraint surface ^ECS), which serves to 
upperbound dependent iterate excursions during successive iterations 
The second corrector phase of the solution strategy lies in the 
use of an energy constraint to scale the generation of successive 
iterates so as to maintain the appropriate form of convergence 
behavior (monotone, oscillating, etc.) associated with the type 
of curvature of the zone of solution space wherein the algorithmic 
tunneling is taking place. The third phase of the solution hierarchy 
involves the use of quality/convergence checks which enable 

various self-adaptive modifications of the algorithmic structure. 

In the sections which follow, detailed discussions are given 
on the classical INR algorithm and its limitations, the develop- 
ment of the various levels of the self adaptive predictor- 
corrector approach as well as the results of several numerical 
examples which demonstrate the capabilities of the new procedure. 
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GOVERNING CLASSICAL INR OPERATOR 
Before overviewing the development of the CINR algorithm, 
it is worthwhile to review the salient features of the INR as well 
as outline several of its more important shortcomings. 


2.1 INR Algorithm 


Assuming that large deformation processes are in effect, 


the virtual work principle takes the following form in Lagrangian 

r2ll 

coordinates namely 


/(6L,.S,j " - 0 


( 2 . 1 ) 


where 6( ), ^ , u^. ^ q^. and R respectively denote the varia- 

r 21 1 

tional operator, 2nd Piola-Kirchhoff stress tensor , 

r 211 

the Lagrangian (Green's) strain tensor , displacement, body 

force and lastly the region occupied by the structure. Introducing 

[151 

the shape function description of displacements , 

U = [N]Y (2.2) 

the following assembled finite element (FE) formulation is 
obtained, that is 

/CB*(Y)]^ Sdv = F(Y). (2.3) 

R ^ ~ 

T 

where ( ) is matrix transposition, and 

[B*] = [B: + CB^(Y)][G] (2.4) 

such that [B], [G] are nonlinear partitions of the strain 

and [N], S and Y respectively represent the shape function, vector 
form of stress tensor and the nodal displacements. 
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Since (2.3) is inherently nonlinear, assuming that the 
material properties can be cast in a tangent stiffness formulation, 
namely 


dS - [Dj][B*]dY 


(2.5) 


then (2.3) can be expanded into a truncated Taylor series to 
yield the following operator: 


AF( Y) ~ [K^(Y)]AY + 0(Y Y) 


( 2 . 6 ) 


Now, expressing (2.6) and (2.3) in algorithmic form yields the 
following I NR operator , that is 


^ I^B*(Yj-'')]‘^ S(Y^''') 


(2.7) 


where i, i, j, []'\ respectively denote the £th 

loadstep, ith iteration, jth intermittent update of the stiff- 
ness, matrix inverse, ith displacement increment of the 2,th loadstep 
and lastly the total nodal displacement and force associated with £th loadstep. 
The convergence criteria typically associated with (2.7) 

involve normed checks of successive displacement increments and 

fg 221 

nodal force imbalances, that is 


I Iay’ 


AY 


i-1 


AY^ 


tol 


(2.8) 




f(y; 






^ < tol 


(2.9) 


where here 


j designates the norm 


n 

Z 

i = l 


hi 


( 2 . 10 ) 
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Most typically, satisfaction of such criteria from increment to 
increment is said to be sufficient to guarantee a convergent 
sol ution . 

To streamline the use of (2.8), the consensus opinion seems 
to advocate that [Kj] be updated and inverted only at the beginning 
of a loadstep^^’ This approach yields the so-called modified 

INR (MINR) operator. As noted earlier, to improve the accuracy/ 
convergence characteristics of such an approach, numerous pseudo 
updates have recently been advocated. Here the BFGS family of 
updates has figured prominently . 

2 . 2 Shortcomi nqs 

While the modified, intermittently /constantly/pseudo (BFGS)^^”^^^ 
updated versions of the INR algorithm converge quadrati cal ly if 
the load increments are sufficiently "small", several shortcomings 
can occur when such is not the case. Additional difficulties are 
also encountered in zones of shallow or changing curvature 
definiteness. This situation can be summarized by the following 
comments: 

i) There is no direct way of preselecting increment size 
as nodal force - deflection space changes curvature; 

ii) There is no direct way of establishing an upper bound 
on the magnitude of the iterated displacement, strain, 
stress. and energy excursions for a given load increment; 

iii) Excessive iterate excursions inherently occur in the 
neighborhood of "shallow" slope zones of the force - 
displacement space, and; 
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iv) Without intermittent or constant updating, the iterated 

version of the MINR can exhibit nonmonotone potentially divergent 

convergence characteristics for monotone increasing/decreasing, 

n 4l 

pos i t i ve/negat i ve definite solution curvatures. 

The excessive iterated dependent variable excursions noted 
above tend to cause drifting from the solution curve. When such 
drift is sufficiently large, depending on the topology of the solu- 
tion space, rather strong nonmonotone type divergence may be ini- 

1141 

tiated as the iteration process continues 
3. CONSTRAINED INR (CINR): PREDICTOR PHASE 

In the context of the remarks made in the previous section, 
it follows that one way to limit the excessive excursions of 
successive iterates is to establish some form of upper bound 
constraint. Riks^^*^^ first considered this approach by developing 
a methodology which features the INR and a special parameter 
controlling the progress of the computation in nodal force-deflection 
space . In geometrical terms, the control parameter selected 
corresponds approximately to the arc length of the equilibrium 
path to be computed. It is introduced into the governing field 
equations. Hence, for a problem with N displacement variables, the 
addition of the contraint equation yields an N + 1 dimensional 
space the solution to which is obtained by the NR method. 

Due to the manner in which Riks^^*^^ casts his constraint 
equation, its direct use with the equilibrium equations tends to 
be somewhat awkward for direct use with the standard FE methodology. 
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To circumvent this difficulty, Crisfield employed the 

n 91 

technique advocated by Batoz and Dhatt *- ■’ for standard displace- 
ment control. Using such an approach, Crisfield recast 

the out of balance force vector as a parametized function of the 
external load vector. Due to the use of an inner product type 
constraint on the allowable displacement iterate excursions, this 
approach enabled the development of an expression which sizes the 
allowable iterative changes in external loading. 

In the subsections to follow, the constrained approach is 
generalized to a more comprehensive and self-adaptive form. This 
will be partly achieved by introducing a more general constraint 
namely the hyper-elliptic constraint surface (HECS). Because of 
the greater adaptability of the HECS, this will enable the CINR to 
act as a refined self-adaptive predictor algorithm. In this context, 
the resulting algorithmic structure will be left flexible enough so 
that in the next section, an energy constraint can be introduced to 
serve in the capacity of the associated corrector algorithm. 

3.1 Hyper- El 1 i psoidal Constraint 
Surface HECS 

As noted earlier, to extend the versatility and adaptability 
of the CINR approach, this paper introduces a more general con- 
straint condition namely the hyper-elliptic constraint surface 
HECSasdefinedbytheexpression 

(3.1) 

where l|*|| designates the Euclidean norm and 
2 
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l|y,li^ = (X (3.2) 

such that referring to Fig. (1), is a warping parameter which 
together with the load increment AF^ defines the curvature/geometry 
of the HECS, while and f^ are respectively the displacement and 
load excursions relative to the starting point of the given load 
increment. Figure 2 schematically illustrates the successive use 
of (3.1) in conjunction with the MINR algorithm. By tying the 
selection to the local curvature of the solution curve, the 

geometry of the HECS can be adaptively updated to improve the 
solution flow. As can be seen from Fig. (2), the HECS itself 
establishes a greatest upper bound possible by the iterative excur- 
sions of the dependent field variables. In particular, for the 
nodal displacements, the maximum allowable excursion for a given 
load increment is defined by the expression 

ik£ii, f (3.3) 

By adjusting AF^ and/or varying bounds can be developed for 
the incremental nodal displacement excursions y^. 

To establish the requisite algorithmic hardware arising from 
the use of the HECS, it follows that outside of turning points and 
bifurcations, there are basically four types of curvature behavior 
associated with the solution curve namely: 

i) Monotone decreasing and positive definite (MDPD); 

ii) Monotone increasing and positive definite (MIPD); 

iii) Monotone decreasing and i ndef i ni te (MDID); and 

iv) Monotone increasing and i ndef i ni te (MIID) 
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Since each places varying demands on the algorithmic apparatus, 
the CINR involving the HECS will be structured to admit all such 
situations. 

A structure generally exhibits MDPD behavior at the outset. 
This case will be used as the starting point of the development. 

Referring to Fig. (2), it follows that using the multidimensional 

starting point of the increment as a local origin of the HECS, 

we have that 

yi = y' - Y° 

ll Iz 

= Yj"' - ll * AYj (3.4) 

t h 

where for the i^ iteration 

= ll - ll~'^ ( 3 . 5 ) 

Similarly, f^ is given by 

ll - <3.6) 

where denotes the incremental loading parameter which is iter- 
atively adjusted until the intersection point of the HECS and the 
solution curve is achieved for the given load increment. 

To start the process, either the MINR, INR, or pseudo INR 
algorithms are used to project the solution curve so as to deter- 
mine its intersection with the HECS. In terms of the modified 
NR strategy defined in Fig. (2), the driving force potential 
enabling this calculation is given by 

Aforce = x] AFj - /{ [B* (vj"' ) ]§ ( Y j’' ) - [B* ( Y J ( Y° ) )d v 

(3.7) 
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Hence, considering the MINR for the moment, 

- [B*(V°)]^S(Y°))dv) 
R 


where [Kj] is updated only at the beginning of the load 
Employing (3.8), (3.4) can be reduced to the form 
T i ■ 1 ■> i u 

where here 


[B*(Y°)fs(vJ))dv 


To obtain the intersection point, substituting (3. 
(3.9) into the relation defining the HECS namely (3.1), 
ing expression is obtained 


Pjdlaj-’ + - 1)(||AF||^>" = 0 

t h 

Solving (3.12) for the incremental loading parameter 
yields 




{ - 


_i-l 

~Z2 


- ‘■'“£2 ^ 


4 


_i-l_i-l 
“^1 “£3 


1/2 

] } 


where here 


_i-l 




(3.8) 

increment . 

(3.9) 

(3.10) 

(3.11) 

) and 

the follow 

(3.12) 

(3.13) 


(3.14) 
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-12 ' *'i (?i5i (2i ) 

(3.15) 

= ^1 (liar’ll,)' - (ii^Eiii,)' 

(3.16) 


The proper sign appearing in (3.13) is chosen by noting that 
for MDPD curvature, the bilinear forms k = 1,2,3 have the 

following types of definiteness namely 


(5 


i-1 
£1 ’ 


_i-l 

-%1 


) > 0 


_i-l 

“J,3 


< 0 


1 , 2 ,.. 


(3.17) 


Here since must itself be positive definite for MDPD solution 
geometries, (3.13) is chosen to take the form 


2E 


J {.-•'■-■I + r("''"'')^ - 

_i-l ^ '-^-£2 ^ ^“jll “£3 -I ^ 


(3.18) 


£1 


In this context, the CMINR is structured as follows 






/{[B*(Yj''')]‘^S(Yj‘‘‘) - CB*(Y°)]’’’s(Yj))dv) (3.19) 

R 


Note for MDPD solution curves, the sequence of successive 
AY^ iterates are themselves positive definite. Contingent on the 
successful satisfaction of the convergence criteria, the global 
external load takes the form 

' fl-l ^ (3.20) 


where I denotes the last iteration count. 
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Because of the foregoing properties, successive iterates 
associated with MDPD portions of the solution curve remain con- 
fined inside the HECS. Such is not the case for MIPD situations. 

As seen in Fig. 3, nonmonotone oscillatory convergence is achieved 
wherein successive iterates alternate between increasingly closer 
inside and outside positions relative to the multidimensional 
intersection of the HECS and the solution curve. 

While the CMINR algorithm defined by (3.19) also applies here, 

ri4l 

since the convergence/quality checks '• •’ used to monitor the state 

of solution development may be keyed in on monotonicity proper- 
ties, it is important to determine the " i n/outs i deness " of succes- 
sive iterates. This enables the determination of a consistent 
convergence process. To check for such properties, the functional 
characteristics of the HECS can be used to establish the in/out- 
sideness of the i^ iterate by evaluations of the following condi- 
tion f 1 ag namely 

- (IIAF^IIJ^ (3.21) 

where 

-i r >0 outside point 

^5, < 0 inside point (3.22) 

Note for such situations, the definiteness characteristics of 
“£k* ^ ~ altered. In particular, since the successive 

solution curvatures are steeper than the initial state, it follows 


that 


„i-l 

“^1 


> 0 ; 


_i-l 


< 0 


(3.23) 


indefinite 


1 = 


1,2,3... 
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FIG.3 Nonmonotone but convergent iterative process associated with 
HECS constrained MINR algorithm in zone of MIPD curvature 
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In the case of MDID local solution behavior, the bilinear 


forms k = 1,2,3 have the following definiteness character- 

istics for Vi namely 


_i-l 

“jn 


> 0 ; 



_i-l 


i ndef i ni te 


1,2,3... 


(3.24) 


Note, in the context of Fig. (4), the force potential driving the 
INR projection of the solution curve into its intersection with the 
HECS is given by the same expression as positive definite situa- 
tions, namely (3.8). Here though, due to the nature of the inter- 
section, the load parameter takes the form 



1 

„i-l 



_i-l 


‘■^-£2 ^ 



_i-l 

“£3 


1/2 

]} 


(3.25) 


Note as with MI P D s i tua t i ons , successive iterates form an oscil- 
latory nonmonotone sequence whose members are alternating inside 
or outside of the HECS. Such properties can be ascertained by 
employing the criterion function defined by (3.21). 

Lastly for MI ID situations described in Fig. (5), all the 
modified algorithms established for the preceeding indefinite 
case also apply here; the only exception being that successive 
iterates display a MDID type behavior and hence remain inside the 
HECS. In this context. 


_i-l 

“£3 


_i-l 

“£2 


) > 0 


< 


0 



1,2, I 


(3.26) 
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FIG. 5 




i>IYI 




IFI, Ill'll 



o IIYI 


«v 


Monotone iterative process associated with HEC3 constrained 
MINR algorithm in zone of MIPD curvature 
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The preceeding algorithms were all developed for some 
general iteration. For the first, several simplifications can 

obviously be made, in particular the load parameter takes the form 





) 


1/2 


) (IIAFJIJ 


(3.27) 


In terms of (3.27), the algorithm defining the successive dis- 
placement iterates for PD and ID situations reduce to the follow- 
i ng form namely 


ay’ = [k^(y;)]~’[ — 


1/2 




(3.28) 





^ (I|AF,||^)^ 


1/2 

(IIAF^II^)^ 


(3.29) 


In the preceeding algorithmic developments, it was tacitly 
assumed that the types of definiteness of the solution curve 
remained unchanged during the successive iterations associated 
with a given load increment. For situations which straddle turn- 
ing points, such is not the case. Since the algorithmic structure 
is different for positive and negative definite situations, some 
provisions must be developed to identify such changes in definite- 
ness so that the proper modifications can be made. To initiate 
adaptive updates of the stiffness as triggered by definiteness 
changes, it is assumed that load incrementing as enforced by the 
HECS is tight enough so that either MDPD or MUD behavior is 
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encountered to the left of turning points. For such situations, 
comparison checks between successive iterates can be used to 
monitor definiteness changes. In this context, accounting for 
the initial curvature of a given load increment, the following 
condition flag can be introduced and monitored namely 

kJ* = sgn(c') sgndiFj-'llj- (3.30) 


where plus to minus sign change can be used to signal updates. 

In the context of (3.30), the algorithm defining takes the form 


1 






•£2 


45 


_i-l _i-l 


H2 “£3 


1/2 

] } (3.31) 


i>l 


An alternative test can be used to trigger the updating of the 
stiffness in the neighborhood of turning points. As seen from Fig. 
(6), successive form a monotone decreasing sequence namely 


X° > X^ > X^ > . . . > X^ > . . . , (3.32) 

While such behavior may initially occur , as seen from 
Fig. (7), passed a certain point, successive x/* values can become 
negative definite namely 

Xj > X^ > . . . > Xj'^ > 0 > Xj > . . (3.33) 

Such a change in definiteness can be used to trigger the update 
process. At such a point, the choice of the proper X^ algorithm 
is keyed in on the definiteness encountered. As an example, for 
turning points which involve transitions from negative to positive 
definite curvature, the monotonicity noted above is reversed. 
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Similar results can also be ascertained by monitoring the 
"above/belowness" relative to the lly|| axis of the HECS, namely 


~ 2 




^ below 


I- 


I F(Y 

' A/ 


i-1 


(3.34) 


; above 


Adjusting for the initial curvature of the given load increment, 
the following condition flag can be used to establish the requi- 
site restructuring of the algorithm, that is 


i-1 


sgn ( l|F(vJ'')ll^ - llf(vj)ll^ ) 

i>l 


(3.25) 


where sign ch^anges signal the need for- stiffness updating such that 

below origin 

(3.36) 


<t> 


1-1 


above origin 

In terms of (3-35), the CMINR algorithm takes the following form: 


= 1 r _-i-l . xi-1 r/-i-lv2 A-i-1 -i-1-]^^^ , 

o-i-1 ^ ~Z2 l(-^2 ' ■ ^-0'’ “'’■5 ^ ^ 


2 - i - 1 


■Z2 ~Z3 


(3.37) 


Note, to keep the above noted algorithmic flow consistent, the + signs appearing 


in (3.27) must be replaced by sgn(C^).This will yield the proper sucession of J 


3. 2 Adaptive Warp of HECS 

To establ is h , the local curvature of the force displace- 
ment space is required. In this context, the curvature parameter 

[91 

of Bergan et al. '■ ■* is particularly useful as it represents a 
measure of the local definiteness (positive or indefinite). For 
the present purposes, to establish such a relation, assuming that 
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AF^ is a 

constant , 

then F^ is defined by the s i ngl e. parameter 

relation 





(3.38) 

where 

^ I 

k = l 


h - 

(3.39) 


In terms of (3.39), the curvature parameter is obtained by taking 
the ratio of the inner products of AF and the derivative of the 
nodal displacement via A|^ evaluated at k = 1 and 1 respectively. 
This yields the expansion 


Cr = (iF)'^ ^ (Y)|Aj.,/((AF)'' ^ (Y)|A,) (3.40) 

[231 

where employing backward finite differences '■ , the foregoing 

derivatives can be approximated by 


dA 








+ o((x‘.,)^) 


(3.41) 


In terms of (3.41), (3.40) reduces to 

-X. (AF)'^AY^_^x{ / ((AF)"^AY^xJ_^ ) (3.42) 


such that AY^ and 4Y^_.| represent the total variations in nodal 
displacements associated with the first and (5,-1)^^ load increments. 

The curvature parameter can be further modified by noting 
that for small enough excursions, it follows that 

xj AF ^ CKT(Yk_i)]AYk (3.43) 


hence 
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" (av,)'^[k^(y?)]ay, xJ_, 


(3.44) 


where here, the denominator is a direct measure of the incremental 
energy stored during the first load step while the numerator 
denotes the second variation of energy associated with the (£,-1)^*^ 
load increment. 

0 

The parameter can be used to scale To start this 

development, it follows that during the initial stages of any 
loading process, only modest changes typically occur in [Kj] hence 
few iterations occur during say the first increment. Thus 

AY^^-[Kt(0)]'^AF (3.45) 

or in a normed sense 

I |AY^ I I 2 < M[Kt(P)]’''aF| I 2 (3.46) 


Recalling the HECS, it follows that the upper bound value of AY.j 
is given by 

1 I Is ; J- 1 1'^?! I 2 

f-^0 ^ 


and hence. 




— AF 


(3.48) 


Comparing (3.46) and (3.47), it follows that a good initial 
value of can be taken as 


p ^ (initial) = 


a 


(3.49) 
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where 

Nj = llAFIIj/ ||[K^(0 )]-'aF||2 (3.50) 

such that a is a user selected parameter which enables an expansion 
or contraction capability for the HECS. Now as we proceed to suc- 
cessive load steps, must be scaled to reflect potential curva- 
ture changes in the force-deflection space. Since = 1 , this can 
be achieved by letting 


/ p ^ \ 3 


(3.51 ) 


where B enables the user to vary the influence of the curvature 
parameter in defining the warping of the HECS. 


4. ENERGY CONSTRAINT: CORRECTOR PHASE 

As noted earlier, for the present purposes the CMINR is 
employed in the manner of a predictor algorithm. To correct 
the results arising from this stage of calculation, a strain energy 
constraint will be employed to enforce the proper type of mono- 
tonicity of successive solution iterates. This is achieved by 
upper bounding the admissible strain energy excursion by scaling 
the variation of load and deflection during the iteration process. 
Such scaling can either be based on worst case individual element 
constraint tests or on an overall global check. If the check is 
failed, to provide for the foregoing scaling, the HECS is shrunken 
so as to maintain the requisite convergence characteristics. 

To initiate the development, a workable expression must be 
obtained for successive strain energy excursions generated during 
the iterative process. In this context, a trapezoidal approximation 
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is employed to evaluate the incremental area "under" the solution 

t h 

curve. Specifically the energy accumulated during the i^ iter- 
ation of the load increment takes the following form namely 

^^1 = 1 ° ((MayJM2)^) (4.1) 

where 

F(Y^'"') = / [B*(Yj^^’'')]'^S(Yj^''') dv (4.2) 

R 

F(Y^) = / [B*(yJ)]’^S(yJ) dv (4.3) 

R 


To achieve the requisite scaling of the governing field vari- 
ables, yJ is recast as follows 



■< A V ■> 

^ It 


(4.4) 


where the scaling parameter must be chosen to enforce the follow- 
ing energy constraint namely 


AeJ < ej^ ; i = 2,3, . . . (4.5) 

such that is a user selected parameter which can either loosen 
or tighten the monotonicity requirements. Hence, once e^^ is 
selected, (4.1) and (4.5) lead to the requisite value of x^- If 
terms of x^ » HECS can be warped in the abscissa dimension by let- 
ting . This effectively reduces its size thereby providing 

a tighter bound on successive AY^. 

To obtain the foregoing scaling, x^ must be extracted from 
(4.1) and (4.5). In this context, since F(Y^) is dependent on the 
disposition of the energy constraint/scaling parameter x^» in terms 
of (4.4), (4.3) can be recast as follows namely 
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E (ll) = 

R 

(xl) / ([S]^[B„(iY^)]^tD^][B*(Yj-'')] + 

R 

[b*(y^-‘‘)]'^[d^][b^(ay;)][g]) dv ay[ + 


i 2 i T .. ^ 

(x^) / [G] [B^(AY^)] [D^][B^(AY})][G] dv AYJ^ (4.6) 

R 

or in approximate form by 

F(Yj) = F(vj"'') + xj 

R 

Employing (4.6), the energy stored during the ith iteration can 
be written in the form 


(xj)^ /([G]'^[B„(iYj)]'^[D^][B*(Yj''')] + 

R 

[B*(Y’’')]’’[D.][B|,(iY’)][G]dvAY* + 

(x{)^/[G]'^[B„(AYi)]’’[DT][BjAY*)][G]dvAY’ 

R 


Rearranging (4.8), we have that 


= 


X ^ + 

s, 1 


(xj) 


^ r' + 

Fjt2 


(x;> 


3 r + 
^£3 


(x^) 


(4.8) 


(4.9) 
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where the various bilinear coefficients take the form 




r’2 = 


-i-K-,1 


R 


(4.10) 

(4.11) 

(4.12) 


■"m = l-iyj /[G]^[B„(iYj)]'^[D^][B„(iyJ)][G]dviYi (4.13) 

R 

i 2 

Truncating (4.8) to 0((X^) ) or less yields the following more 
tractable algorithmic expression for AE^, that is 

0((x[)^) (4.14) 


Now, enforcing the energy constraint defined by (4.5), the 
following general and reduced polynomial expressions are obtained 
for that is 

(x^)^r^4 + (xj)^rj3 + ■ ®R ^ ° (^-■'5) 


or more simply 


(xh^rj2 + xjr' 


£ £1 


i 0 


For simplicity, solving (4.16) for yields 


1 


£ 2p^ 

£2 


-rj, + + 4e„ E 


i-ipi 1-.T 

£ ^£2-' > 


(4.16) 


(4.17) 
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where for PD situations r ^2 ^ noted earlier, defined 

by (4.17) can be used to resize the HECS thereby providing for a 
tighter bound on successive iterations. Having now obtained the 


proper scaling, the energy stored during the Z 
given by 

T 

,1 


th 


■Jltot 


I 

I AE, 
i = l 


load increment is 


(4.18) 


where here AE^ is defined by (4.9) such that local MDPD solution 
curvature is assumed. In such situations, it follows that 

AeJ > 0 for Vi (4.19) 


Similar monotone behavior of the energy increments is also noted 
for MUD solution curvature. 

In the case of MIPD and MDID curvatures, since successive 
iterates form an oscillatory nonmonotone sequence, the energy 
increments themselves give rise to an alternating sequence of 
positive and negative definite terms. For such a situation, the 
specific definiteness of successive energy increments is defined 
as foil ows : 



> 0 if < 0, > 0 

< 0 if > 0, < 0 


5. SUMMARY AND DISCUSSION OF NUMERICAL EXPERIMENTS 


(4.20) 


The overall algorithmic flow associated with the predictor- 
corrector procedure is performed in several main steps. These 
i ncl ude : 

i) The monitoring of the various condition flags; 
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ii) The application of the various predictor-corrector 
constraint algorithms; and lastly, 


iii) The assessment of convergence. 

For the purposes of algorithmic efficiency, the various condi- 
tion flags can themselves be applied in three main levels which 
have the following purposes namely: 

i) To define the geometry of the HECS contingent upon local 
solution curvature; calculate Cj^, y^; 


ii) Locate solution positioning relative to HECS so as to 


enable proper structuring of algorithms; calculate 


„i-l 


$ 


i-l 


and; 

y 


iii) Define conditioning of iterated solution curve via 

several flags noting need for updating and constraint 
tightening; calculate xj etc. 


As noted earlier, depending on the various condition flags, itera- 
tion count and user options, the stiffness may be updated and in- 
verted in the following manner: 

i) Preferential local updates of highly nonlinear elements 

ii) Standard full global update; 

iii) Pseudo updates (BFGS Broyden DFP 

Huang etc.); 


iv) Update only at start of given load increment loop. 

Such actions are preparatory to the application of the various 
predictor/corrector algorithms. The predictor phase consists of 
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projecting the solution curve via the MINR, INR or pseudo INR 
algorithms to determine its intersection with the HECS. The 
corrector phase employs an energy constraint to enforce the proper 
type of convergence. This is achieved by upper bounding the admis- 
sible energy excursion by scaling the variation of load and deflec- 
tion during the iteration process. Such scaling can either be 
based on worst case individual element constraint tests or an over- 
all global check. For the present purposes, the three phases of 

ri4l 

convergence testing discussed by Radovan '■ ■' are advocated here. 

There consist of: 

i) Displacement/force norm checks; 

ii) Quality of convergence tests; and, 

iii) Nonlinearity checks. 

As a demonstration of the approach developed herein, we con- 
sider the following highly nonlinear numerical experiments, namely: 

i) Stretching of a rubber sheet; 

ii) • Large deformation loading of a spherical cap; and, 

iii) Pre- and postbuckling of a centrally loaded arch. 

These problems were chosen to illustrate the predictor-correctors 
capability and efficiency to handle varying types of kinetic, 
kinematic and material induced nonl i neari ty . To enable the calcula- 
tions, special predictor-corrector "plug ins" were developed for 
the ADINA code of Bathe 

To start, the stretching of a rubber sheet is treated first. 
This problem involves both large deformation kinematics and kinetics 
as well as significant material nonlinearity of the Mooney-Ri vl i n 
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type [28], Figure (8) illustrates the geometry, material prop- 
erties as well as the FE mesh used to simulate the problem. Based 
on the use of 2-D plane stress 8 node isoparametric elements. Fig- 
ures (9 and 10) show various aspects of the response behavior of 
the rubber sheet to wide ranging loads. In addition, Figure (9) 
also lists a comparison of the required number of iterations for 
the MINR and predictor-corrector algorithms over the same load range. 
As can be seen, for the given problems, the current approach is more 
efficient. In particular as seen from Figure 9 a 40% improvement is 
achieved for the given problem. This follows from the fact that the 
HECS tends to generate a larger driving force potential over the 
classic INR for the same size load step. Because of this, fewer 
iterations are required. More importantly is the fact that the en- 
tire iteration process is automatic. The only data needed is the 
final load step. Once specified, the load stepping becomes self- 
adaptive. Note, while a, 3 and e^ are user selectable, for all 
the problems considered herein, unity values proved to yield satis- 
factory results. 

Note while the rate of convergence can be modified by changing 
the various conditioning parameters, due to the constraining nature 
of the predictor-corrector algorithm, "unbounded" iterate excursions 
are precluded from occurring. Because of this, unlike the INR algor- 
ithm which yields strongly divergent and unstable successive iterates 
when excessive load steps are employed, the current approach tends to 
yield a stable solution even when a relatively large HECS and loose 
energy constraint are employed. Whatever solu.tion drift that might 
occur is entirely removed by only moderate tightening of the con- 
straints. This strongly stable characteristic makes the predictor- 
corrector algorithm more forgiving as to conditioning choices. 











o 

CD 


( ) ELEMENT NUMBER 



10 Strain energy space of rubber sheet 
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In terms of the spherical cap problem defined in Fig. 11, 

Fig. 12 clearly demonstrates the foregoing behavioral character- 
istics. In particular, as the HECS is tightened, the correct 
limiting behavior is obtained. Note the other results [27] juxta- 
posed on this Figure were obtained through the use of the INR 
wherein iteration was suspended and hence represents essentially 
a straight Euler-Cauchy type incrementation without regard to un- 
balance loads. When iteration is readmitted into the calculations, 

i 

the INR yields highly unstable and divergent solution behavior. 

This is a direct outgrowth of the fact that for the given cap geom- 
etry, while the global load deflection characteristics show posi- 
tive definite behavior, significant unloading occurs locally. As 
seen in Fig. (13), the slopes of the local element energy-load para- 
meter space undergo fluctuations in definiteness. Because of this, 
the overall stiffness can exhibit local "shallowness" hence leading 
to anomalous excursions of the nodal displacements of a given element. 
For the classic INR type operator, such local overshoot tends to grow 
in magnitude as well as spread to neighboring elements ultimately 
leading to a globally divergent solution. For the current approach 
such behavior is completely eliminated by the use of the HECS and 
energy constraint. Because of this, successive iterations can be 
used to eliminate any load imbalances and hence drift. 

Note, for the results depicted in Figure 12, the CINR gener- 
ated results were between 70-80% faster than the standard INR with 
iteration suspended. If small amounts of drift were allowed in the 
CINR, 5% max, the speed of calculation improved to the range 140-100% 
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FIG. 11 FE simulationof spherical cap 
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deflection ratio Wq/H 


FIG. 12 Load deflection curve of spherical cap 
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times faster. If the same was attempted for the INR, as noted 
above, divergent solution behavior was immediately encountered. 

The foregoing speed enhancements are associated with the form of 
nonlinearity treated. Had other types been considered, the speed 
enhancements would have varied depending on the generic curvature 
changes encountered. 

Figure (14) illustrates the geometry and finite element model 
of a centrally loaded shallow arch. The model employs plane stress 
eight node i soparametric elements. As seen from Figure (15), good 
correlation is obtained with previous analytical [29] and experi- 
mental [30] results. The local load/unload characteristics of the 
pre- to postbuckl i ng transitions are clearly seen in Figure (16). 

As with the cap problem, local changes in definiteness occur in the 
energy-load parameter space. For the given arch though, such defin 
iteness fluctuations are significant enough to lead to unloading/ 
reloading in the postbuckl ing zone. 

6. CONCLUSIONS 

In terms of the foregoing, numerical experiments, it follows 
that the predictor-corrector algorithm can handle essentially all 
the types of nonlinearities prevalent to the nonlinear responses 
of structures in a highly efficient and self-adaptive fashion. 

This includes situations which undergo definiteness changes as in 
turning points and bifurcations. Because of the manner of the 
formulation, the procedure is applicable to history dependent sit- 
uations involving creep and plasticity. Lastly, due to the form 




FIG. 14 FE simulation of arch 
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FIG. 15 Load deflection curve of arch 
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CENTRAL DEFLECTION(INCHES) 

FIG. 16 strain energy space of arch 
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of its algorithmic "hardware", it can be easily implanted into 
currently available GP nonlinear codes without any need for major 
architectural modification. 
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Nomencl ature 

Abbreviation 

Meani n 


BFGS 

DFP 

CINR 

CMINR 

FE 

HECS 

GP 

ID 

INR 

MDID 

MDPD 

MUD 

MIPD 

MINR 

PD 


Broyden-Fl etcher-Go 1 dfarb-Shanno update 
Davidon, Fletcher, Powell method 

Constrained Incremental Newton Raphson Method 

Constrained modified Incremental Newton Raphson 
Method 

Finite Element 

Hyper ellipsoidal constraint surface 

General Purpose 

Indefinite 

Incremental Newton Raphson 
Monotone decreasing indefinite 
Monotone decreasing positive definite 
Monotone increasing indefinite 
Monotone increasing positive definite 
Modified Incremental Newton Raphson 
Positive definite 





Nomenclature 


Symbol 


i-1 


£ 

[B] 

[Bn(Y)l [G] 


[B*(Y)1 


[Oj] 


AE, 






Meam' 


Intercept of MINR extrapolation 
solution curve in [ | y| | ,X | | AF | | 


Linearmatrix coefficient of Y defining strain 


Nonlinear matrix coefficient of Y defining 
strain 


Matrix coefficient of Y defining variation 
in strain 


Slope of MINR extrapolation of solution curve 


in I I y| 1 I 1 AF 1 space 


Curvature parameter 
Material stiffness 
Allowable energy ratio 
Energy Increment 
Nodal force 

Increment in nodal force 


[K-r(Y)],[K^] 

L . ,.L 
1 j 

[N] 

f^s 

0 ( ) 

q.,Q 

1 ^ 


^i j 


Load excursions relative to starting point 
of given increment 

Tangent stiffness 

Update parameter 

Lagrangian strain tensor 

Shape function 

Normed quantity use to define 
On the order of ( ) 

Body force 

Initial region occupied by structure 
2nd Piola Kirchoff 


stress tensor 
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